1 Example: Monthly electricity sales for Virginia

1.1 Extract data from remote database

esales <- dbGetQuery(db,'SELECT * from eia_elec_sales_va_all_m') # SQL code to retrieve data from a table in the remote database
# str(esales)
esales <- as_tibble(esales) # Convert dataframe to a 'tibble' for tidyverse work
# str(esales)
# Reference: https://arrow.apache.org/docs/r/
# if(!('arrow' %in% installed.packages())) install.packages('arrow')
library(arrow)
write_feather(esales, "esales.feather")
# Close connection -- this is good practice
dbDisconnect(db)
dbUnloadDriver(db_driver)

1.2 Briefly characterize the dataset

library(arrow)
esales <- read_feather("esales.feather")
print(esales)    # print the data as a table
summary(esales)  # compute basic summary statistics about the data
     value            date                 year          month       
 Min.   : 7153   Min.   :2001-01-01   Min.   :2001   Min.   : 1.000  
 1st Qu.: 8200   1st Qu.:2005-11-01   1st Qu.:2005   1st Qu.: 3.000  
 Median : 9019   Median :2010-09-01   Median :2010   Median : 6.000  
 Mean   : 9093   Mean   :2010-08-31   Mean   :2010   Mean   : 6.425  
 3rd Qu.: 9885   3rd Qu.:2015-07-01   3rd Qu.:2015   3rd Qu.: 9.000  
 Max.   :11724   Max.   :2020-05-01   Max.   :2020   Max.   :12.000  
# References: https://www.tidyverse.org/, https://dplyr.tidyverse.org/

esales %>%
  filter(year == 2019) %>%
  filter(value > 9000) %>%
  print()

esales %>%
  group_by(month) %>%
  summarise(mean = mean(value)) -> mean_esales_by_month

esales %>%
  mutate(sales_TWh = value/1000) %>%
  select(-value)
  
# filter(data object, condition) : syntax for filter() command

1.3 Plot the time series

#Reference: https://ggplot2.tidyverse.org/

ggplot(data=esales, aes(x=date,y=value)) + 
  geom_line() + xlab("Year") + ylab("Virginia monthly total electricity sales (GWh)")

1.4 Convert the data frame into a time series tsibble object

library(lubridate) # Make it easy to deal with dates

esales %>% filter(month==3)                   # These three lines of code
esales %>% filter(month(date)==3)             #   all do
esales %>% filter(lubridate::month(date)==3)  #   the same thing.

# We don't have to keep the 'year' and 'month' column: can recover them if needed

esales %>%
  select(date, sales_GWh = value) -> esales_tbl

print(elsales_tbl)
# install.packages("tsibble")
library(tsibble) # Reference: https://tsibble.tidyverts.org/articles/intro-tsibble.html

esales_tbl %>% as_tsibble(index = date) -> elsales_tbl_ts

print(elsales_tbl_ts)

1.5 Make some plots

1.5.1 Make a histogram of monthly sales

hist(elsales_tbl_ts$sales_GWh, breaks=40)

1.5.2 Make a seasonal plot

elsales_tbl_ts %>% 
  feasts::gg_season(sales_GWh, labels = "both") + ylab("Virginia electricity sales (GWh)")
# install.packages("feasts"), Reference: https://feasts.tidyverts.org/
library(feasts)

elsales_tbl_ts %>% 
  mutate(Month = tsibble::yearmonth(date)) %>% 
  as_tsibble(index = Month) %>%
  select(Month,sales_GWh) -> vaelsales_tbl_ts

print(vaelsales_tbl_ts)
vaelsales_tbl_ts %>% gg_season(sales_GWh, labels = "both") + ylab("Virginia electricity sales (GWh)")

vaelsales_tbl_ts %>% 
  gg_subseries(sales_GWh)

vaelsales_tbl_ts  %>% filter(month(Month) %in% c(3,6,9,12)) %>% gg_lag(sales_GWh, lags = 1:2)


vaelsales_tbl_ts  %>% filter(month(Month) == 1) %>% gg_lag(sales_GWh, lags = 1:2)

vaelsales_tbl_ts %>% ACF(sales_GWh) %>% autoplot()

# if(!('fpp3' %in% installed.packages())) install.packages('fpp3')
library(fpp3)
# decompose(vaelsales_tbl_ts)
vaelsales_tbl_ts %>%
  model(STL(sales_GWh ~ trend(window=21) + season(window='periodic'), robust = TRUE)) %>%
  components() %>%
  autoplot()

vaelsales_tbl_ts %>%
  mutate(ln_sales_GWh = log(sales_GWh)) %>%
  model(STL(ln_sales_GWh ~ trend(window=21) + season(window='periodic'),
    robust = TRUE)) %>%
  components() %>%
  autoplot()

vaelsales_tbl_ts %>%
  features(sales_GWh, feat_stl)
vaelsales_tbl_ts %>%
  features(sales_GWh, feature_set(pkgs=\feasts\))

2 Example: Australian production

# install.packages('tsibbledata')
library(tsibbledata)

aus_production

aus_production %>% gg_season(Electricity)

aus_production %>% gg_season(Beer)

3 Example: Gross Domestic Product data

3.1 Exploratory data analysis

library(tsibbledata) # Data sets package

print(global_economy)
global_economy %>% filter(Country==\Sweden\) %>% print()
# A tsibble: 58 x 9 [1Y]
# Key:       Country [1]
   Country Code   Year          GDP Growth   CPI Imports Exports Population
   <fct>   <fct> <dbl>        <dbl>  <dbl> <dbl>   <dbl>   <dbl>      <dbl>
 1 Sweden  SWE    1960 14842870293.  NA     9.21    23.4    23.0    7484656
 2 Sweden  SWE    1961 16147160123.   5.68  9.41    21.7    22.3    7519998
 3 Sweden  SWE    1962 17511477311.   4.26  9.86    21.4    21.9    7561588
 4 Sweden  SWE    1963 18954132366.   5.33 10.1     21.5    21.9    7604328
 5 Sweden  SWE    1964 21137242561.   6.82 10.5     21.9    22.3    7661354
 6 Sweden  SWE    1965 23260320646.   3.82 11.0     22.5    21.9    7733853
 7 Sweden  SWE    1966 25302033132.   2.09 11.7     21.9    21.4    7807797
 8 Sweden  SWE    1967 27463409202.   3.37 12.2     21.0    21.1    7867931
 9 Sweden  SWE    1968 29143383491.   3.64 12.5     21.6    21.6    7912273
10 Sweden  SWE    1969 31649203886.   5.01 12.8     23.0    22.8    7968072
# … with 48 more rows
global_economy %>%
  filter(Country==\Sweden\) %>%
  autoplot(GDP) +
  ggtitle(\GDP for Sweden\) + ylab(\$US billions\)

3.2 Fitting data to simple models

global_economy %>% model(trend_model = TSLM(GDP ~ trend())) -> fit
fit
fit %>% filter(Country == \Sweden\) %>% residuals()
fit %>% filter(Country == \Sweden\) %>% residuals() %>% autoplot(.resid)

3.2.1 Work with ln(GDP)

global_economy %>%
  filter(Country==\Sweden\) %>%
  autoplot(log(GDP)) +
  ggtitle(\ln(GDP) for Sweden\) + ylab(\$US billions\)

global_economy %>%
  model(trend_model = TSLM(log(GDP) ~ trend())) -> logfit
logfit %>% filter(Country == \Sweden\) %>% residuals() %>% autoplot()

global_economy %>% model(trend_model = TSLM(log(GDP) ~ log(Population))) -> fit3
fit3 %>% filter(Country == \Sweden\) %>% residuals() %>% autoplot()

4 Producing forecasts

fit %>% forecast(h = \3 years\) -> fcast3yrs

fcast3yrs
fcast3yrs %>% filter(Country == \Sweden\, Year == 2020) %>% str()
fable [1 × 5] (S3: fbl_ts/tbl_ts/tbl_df/tbl/data.frame)
 $ Country: Factor w/ 263 levels \Afghanistan\,..: 232
 $ .model : chr \trend_model\
 $ Year   : num 2020
 $ GDP    : dist [1:1] 
  ..$ 3:List of 2
  .. ..$ mu   : num 5.45e+11
  .. ..$ sigma: num 5.34e+10
  .. ..- attr(*, \class\)= chr [1:2] \dist_normal\ \dist_default\
  ..@ vars: chr \GDP\
 $ .mean  : num 5.45e+11
 - attr(*, \key\)= tibble [1 × 3] (S3: tbl_df/tbl/data.frame)
  ..$ Country: Factor w/ 263 levels \Afghanistan\,..: 232
  ..$ .model : chr \trend_model\
  ..$ .rows  : list<int> [1:1] 
  .. ..$ : int 1
  .. ..@ ptype: int(0) 
  ..- attr(*, \.drop\)= logi TRUE
 - attr(*, \index\)= chr \Year\
  ..- attr(*, \ordered\)= logi TRUE
 - attr(*, \index2\)= chr \Year\
 - attr(*, \interval\)= interval [1:1] 1Y
  ..@ .regular: logi TRUE
 - attr(*, \response\)= chr \GDP\
 - attr(*, \dist\)= chr \GDP\
 - attr(*, \model_cn\)= chr \.model\
fcast3yrs %>% 
  filter(Country==\Sweden\) %>%
  autoplot(global_economy) +
  ggtitle(\GDP for Sweden\) + ylab(\$US billions\)

4.1 Model residuals vs. forecast errors

Model residuals:

Your data: \(y_1, y_2, \ldots, y_T\)

Fitted values: \(\hat{y}_1, \hat{y}_2, \ldots, \hat{y}_T\)

Model residuals: \(e_t = y_t - \hat{y}_t\)

Forecast errors:

augment(fit)
augment(fit) %>% filter(Country == \Sweden\) %>%
  ggplot(aes(x = .resid)) +
  geom_histogram(bins = 20) +
  ggtitle(\Histogram of residuals\)

4.2 Are the model residuals auto-correlated?

augment(fit) %>% filter(Country == \Sweden\) -> augSweden

augSweden %>%
  ACF(.resid) %>%
  autoplot() + ggtitle(\ACF of residuals\)

augment(fit3) %>% filter(Country == \Sweden\) -> augSweden3

augSweden3 %>%
  ACF(.resid) %>%
  autoplot() + ggtitle(\ACF of residuals\)

5 Example: GDP, several countries

library(tsibbledata) # Data sets package

nordic <- c(\Sweden\, \Denmark\, \Norway\, \Finland\)

(global_economy %>% filter(Country %in% nordic) -> nordic_economy)
nordic_economy %>% autoplot(GDP)

fitnord <- nordic_economy %>%
  model(
    trend_model = TSLM(GDP ~ trend()),
    trend_model_ln = TSLM(log(GDP) ~ trend()),
    ets = ETS(GDP ~ trend(\A\)),
    arima = ARIMA(GDP)
  )

fitnord
fitnord %>%
  select(arima) %>%
  coef()

Denmark: ARMA(1,1)

Finland: MA(2)

Norway: MA(1)

Sweden: MA(2)

nordic_economy %>%
  model(arima_constrained = ARIMA(GDP ~ pdq(1,0,2))) %>% select(arima_constrained) %>% coef()
fitnord %>% coef() 
fitnord %>%  glance()  
fitnord %>% filter(Country == \Denmark\) %>% select(arima) %>% report()
Series: GDP 
Model: ARIMA(1,1,1) 

Coefficients:
          ar1     ma1
      -0.3898  0.7240
s.e.   0.2061  0.1434

sigma^2 estimated as 2.407e+20:  log likelihood=-1417.5
AIC=2840.99   AICc=2841.45   BIC=2847.12
fitnord %>%
  accuracy() %>%
  arrange(Country, MPE)
# ETS forecasts
USAccDeaths %>%
  ets() %>%
  forecast() %>%
  autoplot()
str(taylor)
plot(taylor)

6 References

LS0tCnRpdGxlOiAgICAgIkV4cGxvcmF0b3J5IGFuYWx5c2lzIG9mIHRpbWUgc2VyaWVzIGRhdGE6IEV4YW1wbGVzIgppbnN0aXR1dGU6ICJTWVMgNTU4MSBUaW1lIFNlcmllcyAmIEZvcmVjYXN0aW5nLCBTcHJpbmcgMjAyMSIgCmF1dGhvcjogICAgICJJbnN0cnVjdG9yOiBBcnRodXIgU21hbGwiCmRhdGU6ICAgICAgICJWZXJzaW9uIG9mIGByIFN5cy5EYXRlKClgIgpvdXRwdXQ6IAogIGh0bWxfbm90ZWJvb2s6CiAgICBudW1iZXJfc2VjdGlvbnM6IHRydWUKICAgIHRvYzogeWVzCiAgICB0b2NfZGVwdGg6IDQKICAgIGNvZGVfZm9sZGluZzogc2hvdyAjIG9wdGlvbnM6IHNob3csIGhpZGUKICAgIGZpZ19jYXB0aW9uOiB5ZXMKICAjIGh0bWxfZG9jdW1lbnQ6CiAgIyAgICAgICBrZWVwX21kOiB5ZXMKICAjIHBkZl9kb2N1bWVudDogZGVmYXVsdApiaWJsaW9ncmFwaHk6IC9Vc2Vycy9BcnRodXIvR2l0UmVwb3MvVGVhY2hpbmcvdGltZS1zZXJpZXMvdHNlcmllcy5iaWIKbGluay1jaXRhdGlvbnM6IHllcwotLS0KCmBgYHtyIHNldCB1cCBjb2RpbmcgZW52aXJvbm1lbnQsIGluY2x1ZGU9RkFMU0UsIG1lc3NhZ2U9RkFMU0V9CiMgbGlicmFyeShkcGx5cikgLS0gZG9uJ3QgbmVlZCB0aGlzIGlmIHlvdSBhcmUgbG9hZGluZyB0aGUgZW50aXJlICd0aWR5dmVyc2UnIHN1aXRlCmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KGx1YnJpZGF0ZSkgIyBGb3IgZWFzeSBoYW5kbGluZyBvZiB0aW1lLWluZGV4ZWQgb2JqZWN0cwoKIyBpZighKCdmcHAzJyAlaW4lIGluc3RhbGxlZC5wYWNrYWdlcygpKSkgaW5zdGFsbC5wYWNrYWdlcygnZnBwMycpCmxpYnJhcnkoZnBwMykKYGBgCgoKIyBFeGFtcGxlOiBNb250aGx5IGVsZWN0cmljaXR5IHNhbGVzIGZvciBWaXJnaW5pYQoKIyMgRXh0cmFjdCBkYXRhIGZyb20gcmVtb3RlIGRhdGFiYXNlCgpgYGB7ciBvcGVuIGNvbm5lY3Rpb24gdG8gZGF0YWJhc2UsIGV2YWw9RkFMU0UsIGluY2x1ZGU9RkFMU0V9CiMgT3BlbiBjb25uZWN0aW9uIHRvIGEgcmVtb3RlIGRhdGFiYXNlCiMgTWFrZSBzdXJlIHlvdXIgVlBOIG5ldHdvcmsgY29ubmVjdGlvbiBpcyBhY3RpdmUgaWYgbmVlZGVkIQoKIyBpZighKCdSUG9zdGdyZVNRTCcgJWluJSBpbnN0YWxsZWQucGFja2FnZXMoKSkpIGluc3RhbGwucGFja2FnZXMoJ1JQb3N0Z3JlU1FMJykKbGlicmFyeShSUG9zdGdyZVNRTCkKCiMgIm15X3Bvc3RncmVzX2NyZWRlbnRpYWxzLlIiIGNvbnRhaW5zIHRoZSBsb2ctaW4gaW5mb3JtYXRpb24Kc291cmNlKCIvVXNlcnMvQXJ0aHVyL0dpdFJlcG9zL1RlYWNoaW5nL215X3Bvc3RncmVzX2RiX2NyZWRlbnRpYWxzLlIiKQoKIyBPcGVuIGNvbm5lY3Rpb24KZGJfZHJpdmVyIDwtIGRiRHJpdmVyKCJQb3N0Z3JlU1FMIikKZGIgPC0gZGJDb25uZWN0KGRiX2RyaXZlcix1c2VyPXVzZXIsIHBhc3N3b3JkPXBhc3N3b3JkLGRibmFtZT0icG9zdGdyZXMiLCBob3N0PWhvc3QpCnJtKHBhc3N3b3JkKSAKCiMgY2hlY2sgdGhlIGNvbm5lY3Rpb246IElmIGZ1bmN0aW9uIHJldHVybnMgdmFsdWUgVFJVRSwgdGhlIGNvbm5lY3Rpb24gaXMgd29ya2luZwpkYkV4aXN0c1RhYmxlKGRiLCAibWV0YWRhdGEiKQpgYGAKCmBgYHtyIHJldHJpZXZlIGRhdGEgZnJvbSBkYiwgZXZhbD1GQUxTRSwgbWVzc2FnZT1GQUxTRX0KZXNhbGVzIDwtIGRiR2V0UXVlcnkoZGIsJ1NFTEVDVCAqIGZyb20gZWlhX2VsZWNfc2FsZXNfdmFfYWxsX20nKSAjIFNRTCBjb2RlIHRvIHJldHJpZXZlIGRhdGEgZnJvbSBhIHRhYmxlIGluIHRoZSByZW1vdGUgZGF0YWJhc2UKIyBzdHIoZXNhbGVzKQplc2FsZXMgPC0gYXNfdGliYmxlKGVzYWxlcykgIyBDb252ZXJ0IGRhdGFmcmFtZSB0byBhICd0aWJibGUnIGZvciB0aWR5dmVyc2Ugd29yawojIHN0cihlc2FsZXMpCmBgYAoKYGBge3Igc2F2ZSBkYXRhIGluIEFwYWNoZSBBcnJvdyBmb3JtYXQsIGV2YWw9RkFMU0V9CiMgUmVmZXJlbmNlOiBodHRwczovL2Fycm93LmFwYWNoZS5vcmcvZG9jcy9yLwojIGlmKCEoJ2Fycm93JyAlaW4lIGluc3RhbGxlZC5wYWNrYWdlcygpKSkgaW5zdGFsbC5wYWNrYWdlcygnYXJyb3cnKQpsaWJyYXJ5KGFycm93KQp3cml0ZV9mZWF0aGVyKGVzYWxlcywgImVzYWxlcy5mZWF0aGVyIikKYGBgCgpgYGB7ciBjbG9zZSBkYiBjb25uZWN0aW9uLCBldmFsPUZBTFNFfQojIENsb3NlIGNvbm5lY3Rpb24gLS0gdGhpcyBpcyBnb29kIHByYWN0aWNlCmRiRGlzY29ubmVjdChkYikKZGJVbmxvYWREcml2ZXIoZGJfZHJpdmVyKQpgYGAKCiMjIEJyaWVmbHkgY2hhcmFjdGVyaXplIHRoZSBkYXRhc2V0CgpgYGB7ciByZWFkIGluIGRhdGF9CmxpYnJhcnkoYXJyb3cpCmVzYWxlcyA8LSByZWFkX2ZlYXRoZXIoImVzYWxlcy5mZWF0aGVyIikKYGBgCgpgYGB7ciB9CnByaW50KGVzYWxlcykgICAgIyBwcmludCB0aGUgZGF0YSBhcyBhIHRhYmxlCnN1bW1hcnkoZXNhbGVzKSAgIyBjb21wdXRlIGJhc2ljIHN1bW1hcnkgc3RhdGlzdGljcyBhYm91dCB0aGUgZGF0YQpgYGAKCmBgYHtyIHVzZSB0aWR5dmVyc2Ugc3ludGF4IHRvIHBlcmZvcm0gc29tZSBzaW1wbGUgZGF0YSBtYW5pcHVsYXRpb25zfQojIFJlZmVyZW5jZXM6IGh0dHBzOi8vd3d3LnRpZHl2ZXJzZS5vcmcvLCBodHRwczovL2RwbHlyLnRpZHl2ZXJzZS5vcmcvCgplc2FsZXMgJT4lCiAgZmlsdGVyKHllYXIgPT0gMjAxOSkgJT4lCiAgZmlsdGVyKHZhbHVlID4gOTAwMCkgJT4lCiAgcHJpbnQoKQoKZXNhbGVzICU+JQogIGdyb3VwX2J5KG1vbnRoKSAlPiUKICBzdW1tYXJpc2UobWVhbiA9IG1lYW4odmFsdWUpKSAtPiBtZWFuX2VzYWxlc19ieV9tb250aAoKZXNhbGVzICU+JQogIG11dGF0ZShzYWxlc19UV2ggPSB2YWx1ZS8xMDAwKSAlPiUKICBzZWxlY3QoLXZhbHVlKQogIAojIGZpbHRlcihkYXRhIG9iamVjdCwgY29uZGl0aW9uKSA6IHN5bnRheCBmb3IgZmlsdGVyKCkgY29tbWFuZApgYGAKCiMjIFBsb3QgdGhlIHRpbWUgc2VyaWVzCgpgYGB7ciB1c2UgZ2dwbG90MiB0byBnZW5lcmF0ZSBhIHBsb3R9CiNSZWZlcmVuY2U6IGh0dHBzOi8vZ2dwbG90Mi50aWR5dmVyc2Uub3JnLwoKZ2dwbG90KGRhdGE9ZXNhbGVzLCBhZXMoeD1kYXRlLHk9dmFsdWUpKSArIAogIGdlb21fbGluZSgpICsgeGxhYigiWWVhciIpICsgeWxhYigiVmlyZ2luaWEgbW9udGhseSB0b3RhbCBlbGVjdHJpY2l0eSBzYWxlcyAoR1doKSIpCmBgYAojIyBDb252ZXJ0IHRoZSBkYXRhIGZyYW1lIGludG8gYSB0aW1lIHNlcmllcyBgdHNpYmJsZWAgb2JqZWN0CgpgYGB7cn0KbGlicmFyeShsdWJyaWRhdGUpICMgTWFrZSBpdCBlYXN5IHRvIGRlYWwgd2l0aCBkYXRlcwoKZXNhbGVzICU+JSBmaWx0ZXIobW9udGg9PTMpICAgICAgICAgICAgICAgICAgICMgVGhlc2UgdGhyZWUgbGluZXMgb2YgY29kZQplc2FsZXMgJT4lIGZpbHRlcihtb250aChkYXRlKT09MykgICAgICAgICAgICAgIyAgIGFsbCBkbwplc2FsZXMgJT4lIGZpbHRlcihsdWJyaWRhdGU6Om1vbnRoKGRhdGUpPT0zKSAgIyAgIHRoZSBzYW1lIHRoaW5nLgoKIyBXZSBkb24ndCBoYXZlIHRvIGtlZXAgdGhlICd5ZWFyJyBhbmQgJ21vbnRoJyBjb2x1bW46IGNhbiByZWNvdmVyIHRoZW0gaWYgbmVlZGVkCgplc2FsZXMgJT4lCiAgc2VsZWN0KGRhdGUsIHNhbGVzX0dXaCA9IHZhbHVlKSAtPiBlc2FsZXNfdGJsCgpwcmludChlbHNhbGVzX3RibCkKYGBgCgpgYGB7cn0KIyBpbnN0YWxsLnBhY2thZ2VzKCJ0c2liYmxlIikKbGlicmFyeSh0c2liYmxlKSAjIFJlZmVyZW5jZTogaHR0cHM6Ly90c2liYmxlLnRpZHl2ZXJ0cy5vcmcvYXJ0aWNsZXMvaW50cm8tdHNpYmJsZS5odG1sCgplc2FsZXNfdGJsICU+JSBhc190c2liYmxlKGluZGV4ID0gZGF0ZSkgLT4gZWxzYWxlc190YmxfdHMKCnByaW50KGVsc2FsZXNfdGJsX3RzKQpgYGAKCiMjIE1ha2Ugc29tZSBwbG90cwoKIyMjICBNYWtlIGEgaGlzdG9ncmFtIG9mIG1vbnRobHkgc2FsZXMKCmBgYHtyIG1ha2UgYSBoaXN0b2dyYW0gb2YgdGhlIGRhdGF9Cmhpc3QoZWxzYWxlc190YmxfdHMkc2FsZXNfR1doLCBicmVha3M9NDApCmBgYAoKIyMjICBNYWtlIGEgc2Vhc29uYWwgcGxvdAoKYGBge3IsIGV2YWw9RkFMU0V9CmVsc2FsZXNfdGJsX3RzICU+JSAKICBmZWFzdHM6OmdnX3NlYXNvbihzYWxlc19HV2gsIGxhYmVscyA9ICJib3RoIikgKyB5bGFiKCJWaXJnaW5pYSBlbGVjdHJpY2l0eSBzYWxlcyAoR1doKSIpCmBgYAoKYGBge3J9CiMgaW5zdGFsbC5wYWNrYWdlcygiZmVhc3RzIiksIFJlZmVyZW5jZTogaHR0cHM6Ly9mZWFzdHMudGlkeXZlcnRzLm9yZy8KbGlicmFyeShmZWFzdHMpCgplbHNhbGVzX3RibF90cyAlPiUgCiAgbXV0YXRlKE1vbnRoID0gdHNpYmJsZTo6eWVhcm1vbnRoKGRhdGUpKSAlPiUgCiAgYXNfdHNpYmJsZShpbmRleCA9IE1vbnRoKSAlPiUKICBzZWxlY3QoTW9udGgsc2FsZXNfR1doKSAtPiB2YWVsc2FsZXNfdGJsX3RzCgpwcmludCh2YWVsc2FsZXNfdGJsX3RzKQpgYGAKCmBgYHtyICwgd2FybmluZz1GQUxTRX0KdmFlbHNhbGVzX3RibF90cyAlPiUgZ2dfc2Vhc29uKHNhbGVzX0dXaCwgbGFiZWxzID0gImJvdGgiKSArIHlsYWIoIlZpcmdpbmlhIGVsZWN0cmljaXR5IHNhbGVzIChHV2gpIikKYGBgCgoKYGBge3J9CnZhZWxzYWxlc190YmxfdHMgJT4lIAogIGdnX3N1YnNlcmllcyhzYWxlc19HV2gpCmBgYAoKYGBge3IgcGxvdCBsYWdnZWQgdmFsdWVzfQp2YWVsc2FsZXNfdGJsX3RzICAlPiUgZmlsdGVyKG1vbnRoKE1vbnRoKSAlaW4lIGMoMyw2LDksMTIpKSAlPiUgZ2dfbGFnKHNhbGVzX0dXaCwgbGFncyA9IDE6MikKCnZhZWxzYWxlc190YmxfdHMgICU+JSBmaWx0ZXIobW9udGgoTW9udGgpID09IDEpICU+JSBnZ19sYWcoc2FsZXNfR1doLCBsYWdzID0gMToyKQpgYGAKCmBgYHtyfQp2YWVsc2FsZXNfdGJsX3RzICU+JSBBQ0Yoc2FsZXNfR1doKSAlPiUgYXV0b3Bsb3QoKQpgYGAKCmBgYHtyIHBlcmZvcm0gYXV0b21hdGVkIHRpbWUgc2VyaWVzIGRlY29tcG9zaXRpb259CgoKIyBkZWNvbXBvc2UodmFlbHNhbGVzX3RibF90cykKYGBgCgpgYGB7ciBwZXJmb3JtIGFkZGl0aXZlIFNUTCBkZWNvbXBvc2l0aW9uIG9mIHRoZSBWQSBlbGVjdHJpY2l0eSBzYWxlcyB0aW1lIHNlcmllc30KdmFlbHNhbGVzX3RibF90cyAlPiUKICBtb2RlbChTVEwoc2FsZXNfR1doIH4gdHJlbmQod2luZG93PTIxKSArIHNlYXNvbih3aW5kb3c9J3BlcmlvZGljJyksIHJvYnVzdCA9IFRSVUUpKSAlPiUKICBjb21wb25lbnRzKCkgJT4lCiAgYXV0b3Bsb3QoKQpgYGAKCmBgYHtyIHBlcmZvcm0gbXVsdGlwbGljYXRpdmUgU1RMIGRlY29tcG9zaXRpb24gb2YgdGhlIFZBIGVsZWN0cmljaXR5IHNhbGVzIHRpbWUgc2VyaWVzfQp2YWVsc2FsZXNfdGJsX3RzICU+JQogIG11dGF0ZShsbl9zYWxlc19HV2ggPSBsb2coc2FsZXNfR1doKSkgJT4lCiAgbW9kZWwoU1RMKGxuX3NhbGVzX0dXaCB+IHRyZW5kKHdpbmRvdz0yMSkgKyBzZWFzb24od2luZG93PSdwZXJpb2RpYycpLAogICAgcm9idXN0ID0gVFJVRSkpICU+JQogIGNvbXBvbmVudHMoKSAlPiUKICBhdXRvcGxvdCgpCmBgYAoKYGBge3J9CnZhZWxzYWxlc190YmxfdHMgJT4lCiAgZmVhdHVyZXMoc2FsZXNfR1doLCBmZWF0X3N0bCkKYGBgCgpgYGB7cn0KdmFlbHNhbGVzX3RibF90cyAlPiUKICBmZWF0dXJlcyhzYWxlc19HV2gsIGZlYXR1cmVfc2V0KHBrZ3M9ImZlYXN0cyIpKQpgYGAKCiMgRXhhbXBsZTogQXVzdHJhbGlhbiBwcm9kdWN0aW9uCgpgYGB7ciwgd2FybmluZz1GQUxTRX0KIyBpbnN0YWxsLnBhY2thZ2VzKCd0c2liYmxlZGF0YScpCmxpYnJhcnkodHNpYmJsZWRhdGEpCgphdXNfcHJvZHVjdGlvbgoKYXVzX3Byb2R1Y3Rpb24gJT4lIGdnX3NlYXNvbihFbGVjdHJpY2l0eSkKCmF1c19wcm9kdWN0aW9uICU+JSBnZ19zZWFzb24oQmVlcikKYGBgCgojIEV4YW1wbGU6IEdyb3NzIERvbWVzdGljIFByb2R1Y3QgZGF0YQoKIyMgRXhwbG9yYXRvcnkgZGF0YSBhbmFseXNpcwoKYGBge3J9CmxpYnJhcnkodHNpYmJsZWRhdGEpICMgRGF0YSBzZXRzIHBhY2thZ2UKCnByaW50KGdsb2JhbF9lY29ub215KQpgYGAKCmBgYHtyfQpnbG9iYWxfZWNvbm9teSAlPiUgZmlsdGVyKENvdW50cnk9PSJTd2VkZW4iKSAlPiUgcHJpbnQoKQpgYGAKCmBgYHtyfQpnbG9iYWxfZWNvbm9teSAlPiUKICBmaWx0ZXIoQ291bnRyeT09IlN3ZWRlbiIpICU+JQogIGF1dG9wbG90KEdEUCkgKwogIGdndGl0bGUoIkdEUCBmb3IgU3dlZGVuIikgKyB5bGFiKCIkVVMgYmlsbGlvbnMiKQpgYGAKCiMjIEZpdHRpbmcgZGF0YSB0byBzaW1wbGUgbW9kZWxzCgpgYGB7cn0KZ2xvYmFsX2Vjb25vbXkgJT4lIG1vZGVsKHRyZW5kX21vZGVsID0gVFNMTShHRFAgfiB0cmVuZCgpKSkgLT4gZml0CgpmaXQKCgpgYGAKCmBgYHtyfQpmaXQgJT4lIGZpbHRlcihDb3VudHJ5ID09ICJTd2VkZW4iKSAlPiUgcmVzaWR1YWxzKCkKYGBgCgpgYGB7cn0KCmZpdCAlPiUgZmlsdGVyKENvdW50cnkgPT0gIlN3ZWRlbiIpICU+JSByZXNpZHVhbHMoKSAlPiUgYXV0b3Bsb3QoLnJlc2lkKQpgYGAKCiMjIyBXb3JrIHdpdGggbG4oR0RQKQoKYGBge3J9Cmdsb2JhbF9lY29ub215ICU+JQogIGZpbHRlcihDb3VudHJ5PT0iU3dlZGVuIikgJT4lCiAgYXV0b3Bsb3QobG9nKEdEUCkpICsKICBnZ3RpdGxlKCJsbihHRFApIGZvciBTd2VkZW4iKSArIHlsYWIoIiRVUyBiaWxsaW9ucyIpCmBgYAoKYGBge3J9Cmdsb2JhbF9lY29ub215ICU+JQogIG1vZGVsKHRyZW5kX21vZGVsID0gVFNMTShsb2coR0RQKSB+IHRyZW5kKCkpKSAtPiBsb2dmaXQKYGBgCgpgYGB7cn0KbG9nZml0ICU+JSBmaWx0ZXIoQ291bnRyeSA9PSAiU3dlZGVuIikgJT4lIHJlc2lkdWFscygpICU+JSBhdXRvcGxvdCgpCmBgYAoKYGBge3J9Cmdsb2JhbF9lY29ub215ICU+JSBtb2RlbCh0cmVuZF9tb2RlbCA9IFRTTE0obG9nKEdEUCkgfiBsb2coUG9wdWxhdGlvbikpKSAtPiBmaXQzCgpmaXQzICU+JSBmaWx0ZXIoQ291bnRyeSA9PSAiU3dlZGVuIikgJT4lIHJlc2lkdWFscygpICU+JSBhdXRvcGxvdCgpCgpgYGAKCiMgUHJvZHVjaW5nIGZvcmVjYXN0cwoKYGBge3J9CmZpdCAlPiUgZm9yZWNhc3QoaCA9ICIzIHllYXJzIikgLT4gZmNhc3QzeXJzCgpmY2FzdDN5cnMKCmBgYAoKYGBge3J9CgpmY2FzdDN5cnMgJT4lIGZpbHRlcihDb3VudHJ5ID09ICJTd2VkZW4iLCBZZWFyID09IDIwMjApICU+JSBzdHIoKQpgYGAKCmBgYHtyIHZpc3VhbGl6ZSBmb3JlY2FzdHN9CmZjYXN0M3lycyAlPiUgCiAgZmlsdGVyKENvdW50cnk9PSJTd2VkZW4iKSAlPiUKICBhdXRvcGxvdChnbG9iYWxfZWNvbm9teSkgKwogIGdndGl0bGUoIkdEUCBmb3IgU3dlZGVuIikgKyB5bGFiKCIkVVMgYmlsbGlvbnMiKQpgYGAKCiMjIE1vZGVsIHJlc2lkdWFscyB2cy4gZm9yZWNhc3QgZXJyb3JzCgpNb2RlbCByZXNpZHVhbHM6CgpZb3VyIGRhdGE6ICR5XzEsIHlfMiwgXGxkb3RzLCB5X1QkCgpGaXR0ZWQgdmFsdWVzOiAkXGhhdHt5fV8xLCBcaGF0e3l9XzIsIFxsZG90cywgXGhhdHt5fV9UJAoKTW9kZWwgcmVzaWR1YWxzOiAkZV90ID0geV90IC0gXGhhdHt5fV90JAoKRm9yZWNhc3QgZXJyb3JzOgoKYGBge3J9CmF1Z21lbnQoZml0KQpgYGAKCmBgYHtyfQphdWdtZW50KGZpdCkgJT4lIGZpbHRlcihDb3VudHJ5ID09ICJTd2VkZW4iKSAlPiUKICBnZ3Bsb3QoYWVzKHggPSAucmVzaWQpKSArCiAgZ2VvbV9oaXN0b2dyYW0oYmlucyA9IDIwKSArCiAgZ2d0aXRsZSgiSGlzdG9ncmFtIG9mIHJlc2lkdWFscyIpCmBgYAoKIyMgQXJlIHRoZSBtb2RlbCByZXNpZHVhbHMgYXV0by1jb3JyZWxhdGVkPwoKYGBge3J9CmF1Z21lbnQoZml0KSAlPiUgZmlsdGVyKENvdW50cnkgPT0gIlN3ZWRlbiIpIC0+IGF1Z1N3ZWRlbgoKYXVnU3dlZGVuICU+JQogIEFDRigucmVzaWQpICU+JQogIGF1dG9wbG90KCkgKyBnZ3RpdGxlKCJBQ0Ygb2YgcmVzaWR1YWxzIikKYGBgCgpgYGB7cn0KYXVnbWVudChmaXQzKSAlPiUgZmlsdGVyKENvdW50cnkgPT0gIlN3ZWRlbiIpIC0+IGF1Z1N3ZWRlbjMKCmF1Z1N3ZWRlbjMgJT4lCiAgQUNGKC5yZXNpZCkgJT4lCiAgYXV0b3Bsb3QoKSArIGdndGl0bGUoIkFDRiBvZiByZXNpZHVhbHMiKQpgYGAKCgojIEV4YW1wbGU6IEdEUCwgc2V2ZXJhbCBjb3VudHJpZXMKCgpgYGB7cn0KbGlicmFyeSh0c2liYmxlZGF0YSkgIyBEYXRhIHNldHMgcGFja2FnZQoKbm9yZGljIDwtIGMoIlN3ZWRlbiIsICJEZW5tYXJrIiwgIk5vcndheSIsICJGaW5sYW5kIikKCihnbG9iYWxfZWNvbm9teSAlPiUgZmlsdGVyKENvdW50cnkgJWluJSBub3JkaWMpIC0+IG5vcmRpY19lY29ub215KQoKYGBgCgpgYGB7cn0Kbm9yZGljX2Vjb25vbXkgJT4lIGF1dG9wbG90KEdEUCkKYGBgCgpgYGB7cn0KZml0bm9yZCA8LSBub3JkaWNfZWNvbm9teSAlPiUKICBtb2RlbCgKICAgIHRyZW5kX21vZGVsID0gVFNMTShHRFAgfiB0cmVuZCgpKSwKICAgIHRyZW5kX21vZGVsX2xuID0gVFNMTShsb2coR0RQKSB+IHRyZW5kKCkpLAogICAgZXRzID0gRVRTKEdEUCB+IHRyZW5kKCJBIikpLAogICAgYXJpbWEgPSBBUklNQShHRFApCiAgKQoKZml0bm9yZApgYGAKCmBgYHtyfQpmaXRub3JkICU+JQogIHNlbGVjdChhcmltYSkgJT4lCiAgY29lZigpCmBgYAoKCkRlbm1hcms6IEFSTUEoMSwxKQoKRmlubGFuZDogTUEoMikKCk5vcndheTogTUEoMSkKClN3ZWRlbjogTUEoMikKCmBgYHtyfQpub3JkaWNfZWNvbm9teSAlPiUKICBtb2RlbChhcmltYV9jb25zdHJhaW5lZCA9IEFSSU1BKEdEUCB+IHBkcSgxLDAsMikpKSAlPiUgc2VsZWN0KGFyaW1hX2NvbnN0cmFpbmVkKSAlPiUgY29lZigpCmBgYAoKYGBge3J9CmZpdG5vcmQgJT4lIGNvZWYoKSAKYGBgCgpgYGB7cn0KZml0bm9yZCAlPiUgIGdsYW5jZSgpICAKYGBgCgpgYGB7cn0KZml0bm9yZCAlPiUgZmlsdGVyKENvdW50cnkgPT0gIkRlbm1hcmsiKSAlPiUgc2VsZWN0KGFyaW1hKSAlPiUgcmVwb3J0KCkKYGBgCgpgYGB7cn0KZml0bm9yZCAlPiUKICBhY2N1cmFjeSgpICU+JQogIGFycmFuZ2UoQ291bnRyeSwgTVBFKQpgYGAKCgoKYGBge3IsIGV2YWw9RkFMU0V9CiMgRVRTIGZvcmVjYXN0cwpVU0FjY0RlYXRocyAlPiUKICBldHMoKSAlPiUKICBmb3JlY2FzdCgpICU+JQogIGF1dG9wbG90KCkKYGBgCgpgYGB7ciwgZXZhbD1GQUxTRX0Kc3RyKHRheWxvcikKcGxvdCh0YXlsb3IpCmBgYAoKCgojIFJlZmVyZW5jZXMK